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It is becoming ever clearer that the neutrino signal from the next supernova in our Galaxy can reveal 
missing information about the neutrino as well as allowing us to probe the explosion of the star by decoding 
the temporal and spectral evolution of the flavor composition of the signal. But this information may be 
lost if turbulence in the supernova 'depolarizes' the neutrinos so that the observed flux for each flavor 
is an equal mixture of the initial - unencoded - spectra. Determining if depolarization occurs is one of 
the most pressing issues of this field. The most difficult aspect of studying the effect of turbulence upon 
the neutrinos is the lack of any theoretical models that allow us to understand the results of numerical 
studies. This paper makes the suggestion that Random Matrix Theory (RMT) may shine some light in 
this direction and presents support for this the possibility by comparing the distribution of crossing and 
survival probabilities obtained numerically for some 'test case' calculations with the distributions one 
expects from RMT in the calculable limit of depolarization of N neutrino flavors. 
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Introduction 

The subject of neutrino propagation in core-collapse su- 
pernova has been the focus of intense interest in recent 
years with many new effects discovered. From this ever- 
growing body of results we now can expect that the signal 
from the next Galactic supernova to reveal both unknown 
properties of the neutrinos - such as the ordering of the neu- 
trino masses - as well as details of the explosion that would 
allow us to test the basic paradigm. This information is im- 
printed into the neutrino flavor composition of the flux as 
a function of time and energy due to neutrino collective ef- 
fects and the evolving profile of the star as it explodes. For 
a review the reader may consult flllll]. 

But lurking in the shadows is the flavor transforma- 
tion of the neutrino flavor due to the turbulence in the 
supernova. Turbulence in a medium creates additional 
Mikheev-Smirnov-Wolfenstein (MSW) resonances [3, 4] - 
locations where mixing is strongest - which, together with 
the weaker non-resonant fluctuations, can cause flux depo- 
larization [5-8]. Flux depolarization is fatal for any search 
for information in the supernova neutrino signal because 
spectral features are removed or will not be imprinted. De- 
termining if depolarization indeed occurs is a very high pri- 
ority. But the purpose of this paper is not to reconsider the 
potential of the supernova neutrino signal to reveal prop- 
erties of the supernova and/or the neutrino. Instead we 
grapple with the problem that studying and quantifying the 
effect of turbulence in supernova lacks theoretical models 
that can provide insight into the problem. We make the 
suggestion that Random Matrix Theory (RMT) may be the 
doorway to progress in future studies. 

At heart neutrino propagation through a supernova is 
a scattering problem in which we try to relate observed 
neutrino states to their initial states. The applicability of 
RMT to scatteriiigjjroblems has long been recognized: see, 
for example ll9l 4l2ll . That RMT may be useful for neu- 
trino propagation is suggested by the similarity of the non- 



adiabaticity parameter in Kneller & McLaughlin |fl3ll to a 
Breit-Wigner resonance and recalling Ericson's prediction 
[14] from RMT that of fluctuations in nuclear cross sec- 
tions from intermediate states of compound nuclei. If RMT 
is indeed useful then it must be able to predict the distribu- 
tions of survival and crossing probabilities obtained from 
numerical turbulence studies. This is no trivial challenge 
and deciding if RMT has any shot of success is left to the 
future. The view that RMT may be one path to progress 
is supported by its success in one particular case - the de- 
polarized limit - where one can derive the expected distri- 
bution with relative ease. The plan of this paper is thus to 
present the results of a small numerical study of turbulence 
and neutrinos for some selected test cases, to derive the 
expected distributions in the depolarized limit from an en- 
semble of N-flavor random unitary matrices [15], and then 
make the comparison between them pointing out the sim- 
ilarities but also the differences that future, more sophisti- 
cated, applications of RMT will have to predict. 
Turbulence And Supernovae Neutrinos Test Cases 

The quantities we are interested in calculating are the 
survival and crossing probabilities: i.e. Pii'j — s- fi) = Pij 
since the flux of neutrinos emerging from the supernova is 
given by these quantities multiplied by the appropriate ini- 
tial fluxes at the proto-neutron star. Throughout this paper 
we chose to calculate the 'matter' basis probabilities be- 
cause in this basis there are no confusing adiabatic MSW 
transitions. The reader may find the relationship between 
the matter and flavor bases in Kneller & McLaughlin llTll . 
These probabilities can be calculated from the S matrix 
as Pij = The state of the neutrino after travel- 

ing a distance r is related to the initial state at r = by 
|z^(r)) = 5(r, 0) |i^(0)) where the matrix 5'(r, 0) is the 
solution of the equation 



with initial condition that 5(0, 0) = 1. The Hamiltonian 
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H has two components: the vacuum mass term K and the 
potential V that accounts for the effect of matter. We shall 
not consider the case where the neutrino density is so high 
that collective neutrino effects need to be included. The 
vacuum mass term in the flavor basis, K'^^\ is set by the 
mass square differences and the parameters O12, 9i3, O23 of 
the Maki-Nakagawa-Sakata-Pontecorvo (MNSP) mixing 
matrix. The two Majoranna phases of the MNSP matrix are 
irrelevant and we set its CP phase 6 to zero. For this paper 
we adopt dmf^ = 8 x lO^^eV^, \6mU = 3 x 10^^ eV^ 
sin^ 26*12 = 0.83 and sin^ 26*23 = 1 fill]. All values of the 
unknown mixing angle we shall consider will be 'large' 
i.e. they will be above the Dighe & Smirnov |19] thresh- 
old of sin^(20i3) ~ 10^^. If the sign of dml^ is positive 
then the 'hierarchy' - the ordering of the masses - is known 
as 'normal'; if the sign negative then the hierarchy is 'in- 
verted'. The potential affects only the electron neutrinos 
and antineutrinos and, in the flavor basis, the only non- 
zero element of V^^^ is Vjf = y/2GFY^{r)p{r) /ttin 
where Y^{r) is the electron fraction, and rrijq the nucleon 
mass. For simplicity we shall use l^(r) = 0.5. So in 
order to compute the P^ 's we require a density profile 
p{r) for the neutrinos to pass through. To study super- 
nova turbulence p[r) must possess fluctuations but rather 
than adopt density profiles taken from multi-dimensional 
hydrodynamical simulations - which presumably contain 
turbulence already - our approach instead is to start with 
a profile taken from a spherically symmetric simulation, 
i.e. a one dimensional model, that must be turbulence 
free and then add turbulence. In this way we can tune 
the parameters that describe the turbulence plus the un- 
derlying profile is the same for every instantiation so that 
non-turbulent features (such as the shocks) are always the 
same. The density profile we use is the t = 4.5 s snapshot 
of the Q = 3.36 x 10^^ ID hydro described in Kneller, 
McLaughlin & Brockman ifl^ . This particular simulation 
was chosen because it most closely resemble the profiles 
of the 2D model described there. We do not distribute the 
turbulence throughout the entire profile because the profile 
is dissected into distinct regions by the forward and reverse 
shocks. Within each region one would expect the turbu- 
lence, if any, to be quite different. For this paper we shall 
only insert turbulence between the reverse shock at and 
the forward shock at which is the region with largest 
amplitude fluctuations as indicated by multi-dimensional 
supernova simulations ll2l1 - t^ll . The density profile p{r) is 
thus constructed as p{r) = {p{r)) + 5p{r) where {p{r)) 
is the adopted one-dimensional profile from the simulation. 
As is common, we model 5p{r) as 6p{r) = F{r) {p{r)) 
with F{r) a random field with vanishing expectation value. 
For this paper, F{r) is proportional to: a parameter 
which sets the amplitude of the fluctuations, the factor 
tanh((r — rj.)/\) tanh((rs — r)/\) whose purpose is 
to avoid introducing additional discontinuities at and 
Ts and where A is a scale set to A = 100 km, and to 
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FIG. 1 : An instance of the potential (top panel) and the survival 
probability Pn (bottom panel) for a i5 = 25 MeV antineutrino 
as functions of the radius r. The turbulence region is identified by 
the vertical dashed lines in both panels, the H resonance density is 
the horizontal line in the upper panel. The underlying base profile 
from the simulation is shown as the dashed line, the evolution 
through the base profile is the dashed line in the lower panel. The 
fluctuation amplitude is = 0.02, the hierarchy inverted and 
sm2(2e'i3) = 4 X lO"''. 

a Gaussian random field described by a power spectrum 
E[k) given by 

E{k) = {a-l){^y . (2) 

The scale /c^, is set to 7r/(rs — r^) i.e. a wavelength twice 
the distance between the shocks, and we shall use a Kol- 
mogorov spectrum where a = 5/3. The 'Randomiza- 
tion' method found in Majda & Kramer |l20l], using 100 
wavenumbers with a cutoff at fc^, is used to instantiate the 
Gaussian random field. 

Now that we have our density profiles we can calcu- 
late how a neutrino propagates through it for some fixed 
set of oscillation parameters and neutrino energy. The 
method used to solve equation ([T]) is described in Kneller 
& McLaughlin [17]. One finds that the probabilities Pij 
for the most past are constant as a function of the distance 
r and change only in the vicinity of the 'resonances' which 
are the locations where the separations of the eigenvalues 
of H are minimal. Resonances occur at two different den- 
sities: the resonance at high density is known as the H res- 
onance, the one at low density is the L resonance. The L 
resonance always occurs between states Ui and U2 because 
the sign of 5m\2 is known. For a normal hierarchy the H 
resonance occurs between states z/2 and z/3 : for an inverted 
hierarchy it is antineutrino states states Ui and that mix. 
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FIG. 2: A normalized frequency histogram of 1000 calculations 
of Pii for a i5 = 25 MeV antineutrino with the underlying profile 
shown in figure (O, an inverted hierarchy and sin^(2Si3) = 4 x 
10-^ 



^0.5 



Numerical Results 
Inverted Hierarchy 

With the setup complete we present numerical results; we 
consider first the case of an inverted hierarchy. In fig- 
ure ([ll we show the evolution of the survival probability 
Pii = P{vi — 7- i^i) as a function of the radius r through 
one instantiation of the potential. Without fluctuations the 
neutrino experiences three H resonances; when turbulence 
is added the number of MSW resonances increases dramat- 
icafly. Since the value is 'large' we observe that the 
semi-adiabaticity of the additional MSW resonances kicks 
the state of the neutrino at it passes through them. But the 
effect of turbulence is not just the addition of new H res- 
onances because we also observe that non-resonant fluctu- 
ations also cause the probability Pn to evolve. Once the 
neutrino has passed through the entire turbulence region 
we find there is a significant difference between the final 
states of the neutrino through the profiles with and without 
turbulence. This difference is not fixed: if we use a new in- 
stantiation of F we obtain a different result. By repeating 
the calculation with 1000 different instantiations we gener- 
ate figure (O which is a frequency distribution of the final 
state i.e. the state at the edge of the star. The figure in- 
dicates that the distribution of the final states is consistent 
with uniform and the final state is completely uncorrelated 
with the initial state. So for this energy, ^13 and turbulence 
amplitude we have a case where two flavor depolarization 
occurs. 

Normal Hierarchy 

When we switch to a normal hierarchy the H reso- 
nance now affects the neutrinos. For modest turbulence, 
Ct < 0.01, we obtain uniform distributions of final states 
for V2 and v^. But as increases we find something 
new. The reason is that for the normal hierarchy the L 
and H resonances both occur for the neutrinos so large am- 
plitude fluctuations break HL factorization i.e. the L reso- 
nances no longer occur after all the H resonances but rather 
the two are mixed together. An example of a case with 
broken HL factorization can be seen in figure (O. For a 
normal hierarchy the H resonance mixes states 1^2 and 1/3 
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FIG. 3: An example of broken HL factorization. The potential 
(top panel) and neutrino survival probability Pn of a 25 MeV 
neutrino (bootom panel) as a function of distance r using the 
same average profile shown in figure ([TJ. The fluctuation am- 
plitude is = 0.1, the hierarchy normal and S13 = 9°. 

and state Vi is unaffected but for this particular calculation 
the density drops very close to the L resonance at around 
r ~ 50, 000 km whereupon the matter state Ui mixes with 
matter state V2- The density thereafter retums to the H res- 
onance where further mixing between states U2 and oc- 
curs. This mixing between states Vi and 1^2 means that the 
distributions for the final state probabilities are now quite 
different. The distributions of the final states for the neu- 
trinos from 1563 calculations for this hierarchy, = 0.3 
and ii^ = 60 MeV are shown in figure (|4|. Without turbu- 
lence Pii is always unity for the mixing parameters used 
and one notices that this is still a common result. But we 
also see that there appears to be a triangular component in 
both the distributions of Pn and P22 '■ if state i^^ has a trian- 
gular component it is not apparent. We thus complete our 
rudimentary study of turbulence effects upon supernovae 
neutrinos; in the next section we derive the expected dis- 
tributions in the limit of N flavor depolarization and then 
compare them with these numerical results. 
Predictions For Random Unitary S Matrices 

We make the conjecture that in the depolarized limit 
the S{r, 0) matrices have a block structure such that each 
block is a random N-flavor unitary matrix ifisll . From this 
ansatz we can derive the distribution of the probabilities 
Pij = I Sij I ^ in this limit. The N real components, , plus 
the N imaginary components, y^j, of the elements in a row 
or column of the block form a 2N Euclidean space. The 
requirement of unitarity defines a unit sphere in this space 
and since these 2N quantities are identically distributed the 
probability of a particular set of the elements must be uni- 
form over the surface of the sphere. That is, the probability 
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Discussion 

We find that the expected distributions in the depolar- 
ized limit possess some resemblance with the numerical 
results. For the test case results shown in figure (O we have 
achieved two flavor depolarization since the distribution is 
uniform; the results shown in figure appear to be a tran- 
sitional stage between the N=2 and N=3 depolarized limits. 
This resemblance is encouraging and supports our sugges- 
tion that better RMT calculations - obviously not based on 
Dyson's ensemble - may able to predict the distributions 
for all cases. This suggestion is the principal message of 
this paper though we are fully aware that deriving the dis- 
tributions for general cases will certainly be a challenge. 
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FIG. 4: Normalized frequency distributions of the probabilities 
Pii, P22 and P33 of 1563 calculations using the same average 
profile shown in figure IT}. The hierarchy is normal, = 0.3 
and £ = 60 MeV. 

that we are located at {xij,yij,X2j,y2j, • • •} where, say 
for a column Sij = Xij + lyij, 82^ = X2j + iy2j, • • •, is 
simply proportional to the area element dA 

P{xij,yij,X2j,y2], ■ ■ .)d^xd^y oc dA (3) 

(N N \ N N 

1 - E ^1 - E Vl n n dy^, ■ 
1=1 i=l / i=l i=l 

We change variables so that the N independent pairs 

x^j , yi-j are expressed as 

~ \/ Pij cosOij^ yij = \fP\j sin^ij, (5) 
X2j = ^/pTj cos 02j, y2j = \/prj sin 92j, (6) 

and so on. After inserting these new variables into equation 
and integrating over the angles we find that 

i=l / i=l 

(7) 

Thus, the set {Pij, ■ ■ ■ Pn]} are uniformly distributed on 
the surface of a standard — 1 simplex. By integrating 
equation (|7]l over — 1 of the P's and normalizing one 
derives our final result that element Pij must be distributed 
according to 

P{P,,) = {N-l){l-P,,r-' (8) 

For N = 2 the distribution is uniform with mean 1/2 and 
variance 1/12; for = 3 the distribution is triangular 
with variance 1/18. 
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